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Abstract 

We show that time-correlation functions of arbitrary order for any random variable in 
a statistical dynamical system can be calculated as higher-order response functions of the 
mean history of the variable. The response is to a "control term" added as a modification 
to the master equation for statistical distributions. The proof of the relations is based upon 
a variational characterization of the generating functional of the time-correlations. The 
same fluctuation-response relations are preserved within moment-closures for the statistical 
dynamical system, when these are constructed via the variational Rayleigh-Ritz procedure. 
For the 2-time correlations of the moment-variables themselves, the fluctuation-response 
relation is equivalent to an "Onsager regression hypothesis" for the small fluctuations. For 
correlations of higher-order, there is a new effect in addition to such linear propagation of 
fluctuations present instantaneously: the dynamical generation of correlations by nonlinear 
interaction of fluctuations. In general, we discuss some physical and mathematical aspects of 
the Ansatze required for an accurate calculation of the time correlations. We also comment 
briefly upon the computational use of these relations, which is well-suited for automatic 
differentiation tools. An example will be given of a simple closure for turbulent energy 
decay, which illustrates the numerical application of the relations. 
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1 Introduction 



It is well-known that, in statistical equilibrium systems, there are very useful relations be- 
tween two-time correlation functions and mean response functions |l|, Q]. The most well-known 
form of this relation gives the 2-time correlation function in terms of a response function of 
the solution of the microscopic equation of motion to an imposed infinitesimal perturbation, 
when the response is averaged over the equilibrium ensemble. These relations are often called 
"fluctuation-dissipation relations" but we prefer the term fluctuation-response relation (FRR) 
as being more descriptive. A similiar relation has been shown to hold arbitrarily far from 
thermodynamic equilibrium in stochastic dynamical systems described by nonlinear Langevin 
equations Q. In this case, however, the response is to a forcing term added into the Fokker- 
Planck equation rather than to the dynamical equation for individual realizations. The validity 
of this form of the theorem depends upon a correct coupling of the force, which, unfortu- 
nately, requires a knowledge of the steady-state invariant measure. This latter fact makes the 
generalized theorem quite difficult to apply in practice. 

It is the purpose of this work to prove a far-reaching generalization of the fluctuation- 
response relation. Our version of the theorem holds for any (time-dependent) Markov process 
described by a master equation for the distribution function in phase space: 

dtV(x,t) =L(t)V(x,t), (1.1) 

We include in our discussion the limiting case of the Liouville equation for a deterministic 
dynamical system. Our theorem is more similiar to that in |3], since it considers the response 



to a driving or "control" term added into the master equation (1.1) rather than to the equation 
for individual realizations. However, in contrast to that result, the coupling of our control 
field does not require any knowledge of the steady-state measure and is quite easy to write 
down explicitly. Most importantly, all multi-time correlations of any finite order are obtained 
as higher-order response functions to the same control field. Furthermore, the statistics of the 
system need not be those of thermal equilibrium or even stationary in time. The proof of 



the relations is based upon a variational characterization of the generating functional for the 
time-correlation functions, which was established in previous work Here we shall give a 
reformulation of that result which is of interest in its own right, as it considerably simplifies 
and streamlines the analysis in the old work. 

The FRR we derive is, however, prohibitively difficult to apply when ( |1,1[ ) describes a 
spatially-extended system with many degrees-of-freedom. In such cases the master equation 
is a PDE in a huge number of variables, far too many to permit a direct numerical solution. 
The practical use of the FRR in this context will depend upon the employment of moment- 
closure approximations. As we shall show, the FRR remains valid within the moment closures 
when these are formulated variationally via the Rayleigh-Ritz method proposed in M. We shall 
review here the Rayleigh-Ritz approximation, providing some new derivations of the old results 
in addition to establishing the FRR's for the moment closures. 

The plan of this paper is as follows: in Section 2 we discuss the variational approach to 
statistical dynamics. The treatment here will be new in several points and provide significant 
simplifications of that in We then employ the variational apparatus to establish the FRR's. 
In Section 3 we review the Rayleigh-Ritz formulation of moment closure and the FRR's in 
that approximation. We also discuss there the physical significance of the FRR's, relating 
that for the 2-time correlation to a linear "regression hypothesis". The closure FRR for the 
(n + l)-time correlations, with n > 1, contains also an effect from the nonlinear terms in the 
closure, namely, the creation of correlations by interaction of fluctuations. Some physical and 
mathematical properties required of the PDF models employed for moment closure will finally 
be discussed. Section 4 concerns numerical aspects, in particular computationally efficient and 
accurate methods for computing the derivatives required in the FRR's. A simple example of a 
turbulence closure will be used to illustrate the numerical issues. The last Section 5 will contain 
our conclusions from this work. An Appendix is also included which summarizes the results 
of the closure FRR diagrammatically in terms of Feynman-type graphs with propagators and 
vertices generated from the closure. 

3 



2 Variational Formulation &; Fluctuation-Response Relations 

(2.1) Variational Approach to Statistical Dynamics 

Suppose that X(i) is (vector-valued) Markov process, whose distribution P(x, t) at time t is 
governed by the forward Kolmogorov equation or master equation 

dtP{x,t) =L(t)7>(x,t), (2.1) 

with L(t) the instantaneous Markov generator. The random process governed by a stochas- 
tic differential equation is a particular example, for which the generator is the Fokker-Planck 
operator. This includes the degenerate case of a deterministic dynamics, for which the gen- 
erator is the Liouville operator. Observables, or random variables, -4(x, t) evolve under the 
corresponding backward Kolmogorov equation 

d t A(x, t) = -L*{t)A(x, t), (2.2) 

in which L*(t) is the adjoint operator of L(t) with respect to the canonical bilinear form on 
L°° x L 1 , i.e. < A, V >:= J dx A(x)V(x). The backward and forward Kolmogorov equations 
may be simultaneously obtained as Euler-Lagrange equations for stationarity of the action 
functional 

F[A,V] := f h dt < A(t), (d t - L(t))V{t) > (2.3) 
Jti 

when varied over V(t) £ L 1 with initial condition V{ti) = Vq and A(t) € L°° with final condition 
A(t f ) = 1. For details, see (§. 

Let Z(i) := Z(X.(t),t) be a random variable for the system given by the continuous function 
3(x,t). Then, the cumulant generating functional Wz\)n\ is defined as 

W z [h] = log(exp ^ dt h T (t)Z(t))>. (2.4) 

The reth-order multi-time cumulants of Z(t) are obtained from W^[h] by functional differenti- 
ation with respect to the "test history" h(i): 

S n W z [h] 



Cii— i n (^1 1 ••■> ^n) 



Shi^ti) ■ ■ ■ 5h in (t n ) 



(2.5) 
h=0 



It is not hard to check from its definition (O) that W z [h] is a convex functional of h. The 



Legendre dual of this functional is defined to be the effective action of Z(i): 

T z [z] = sup{< h,z > -W z [h]}, (2.6) 

h 

with < h,z >:= / dt h T (t)z(t). It is a generating functional of so-called irreducible correlation 
functions of Z(t): 

(2.7) 



r it t * 6nTz[z] 



5z h (ti) ■ ■ ■ Sz in (t n ) z= -' 

The functional derivatives here are evalutated at the mean history z(i) := (Z(i)). It is not hard 
to check from the definition fl2.6j ) that Tz[z] is a convex, nonnegative functional with a unique 
global minimum (equal to zero) at the mean history z = z. 

There is a useful characterization of the effective action Tz[z] by means of a constrained 
variation of the action T[A, V], which was established in |J]. In fact, 

F z [z] = st.pt. A>P T[A,V] (2.8) 

when varied over the same classes as above, but subject to constraints of fixed overlap 

< A(t),V(t) >= I (2.9) 

and fixed expectation 

< A(t),Z(t)V(t) >= z(t) (2.10) 

for all t € [U,tf]. Note that Z(t) is used to denote the operator (in both L 1 and L°° ) of 
multiplication by 2(x,t). The Euler-Lagrange equations for this constrained variation may be 
obtained by incorporating the expectation constraint ( 2,1C| ) with a Lagrange multiplier h(i). 



The overlap constraint may also be imposed with a Lagrange multiplier X(t), as it was in 
However, it turns out to be advantageous to impose ( |2.9| ) through the definitions 



A(t) := l + [B(t)-(B(t)) t ] 

:= 1+C(t), (2.11) 



with the final conditions B(t f ) = C(t f ) = 0. Note that (B(t)) t ■=< B(t),V(t) > is the ex- 
pectation with respect to the distribution V(t). Hence, the overlap constraint ( |2.9D is satisfied 
when B(t) is varied independently of V(t). The variable C(t) is no longer independent of V(t), 
but must satisfy the orthogonality condition < C(t),V(t) >= 0. The expectation constraint 
must still be implemented by the Lagrange multiplier h(t). In terms of B(t) or C(t) the latter 
constraint is 

z(t) = (Z(t)) t + [(Z(t)B(t)) t -(Z(t)) t (B(t)) t } 

= (Z(t)) t + (Z(t)C(t)) t . (2.12) 

The Euler-Lagrange equations are obtained by varying the action T[A, V] over B(t),V(t) with 
A(t) = 1 + [B(t) — (B(t))t\, incorporating the constraint ( 2.12j ) with the multiplier h(t). A 
straightforward calculation gives 

d t V{t) = L(t)V(t) + h T (t)[Z(t) - (Z(t)) t ]V(t) (2.13) 

and 

d t B(t) + L*(t)B(t) + h T (t)[Z(t)B(t) - (Z(t)) t B(t) - (B(t)) t Z(t)} + h T (t)Z(t) = 0. (2.14) 
Let us introduce the new operator 

L h (t) := L(t) + h T (t)[Z(t) - (Z(t)) t }. (2.15) 
The variational equations are written in terms of this operator as 

dtP(t) = L h (t)V(t) (2.16) 

and 

d t B(t) + L* h (t)B(t) + h T (t)Z(t)[l - (B(t)) t ] = 0. (2.17) 

Using the resulting identity 4r(B(t))t = — h T (t)[l — (B(t)) t ](Z(t)) t , the latter equation can be 
rewritten in terms of C(t) = B(t) — (B(t))t as 

8 t C(t) + L* h (t)C(t) + h T (t)[Z(t) - (Z(t)) t ] = 0. (2.18) 



The action functional may be expressed in terms of C, V as 

T[C,V]=j' f dt <C(t),(d t -L(t))V(t)>, (2.19) 

using L*(t)l = 0. The effective action r#[z] is then obtained by substituting the solutions of 
( 2.1G| ),( p,18| ), when the "control field" h(i) is chosen so that ( 2.1 2| ) reproduces the considered 



history z(i). The quantity z(t) can be seen to be "controllable" by h(i) from Legendre duality 
That is, the control h[t; z] for a specified z(i) is obtained from the minimization of the convex 
function W^[h]— < h, z >. Compare equation ( |2,6| ). Gathering together all of our previous 
discussion we may state: 

Proposition: The effective action of the variable Z(i) is obtained as 



T z [z] = I* dt < C(t), (d t - L(t))V{t) >, (2.20) 
Ju 

where C, V satisfy 

d t V(t) = L h (t)V(t) (2.21) 

and 

d t C(t) + L* h (t)C(t) + h T (t)[Z(t) - (Z(t)) t ] = (2.22) 
with initial-final conditions 

V(ti)=V , C(t/) = 0, (2.23) 
and the value of the control field h is selected to give for all t € [ti,tf] 

(Z(t)) t + (Z(t)C(t)) t = z(t). (2.24) 

(2.2) Fluctuation-Response Relations 

It is not accidental that the same notation h(t) was chosen above for the control field as for the 
argument of the cumulant-generating functional W^[h]. In fact, we shall prove that 

W z [h] = ( S dt h T (t)(Z(t)) t , (2.25) 
Ju 
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using just the solution V(t; h) of the forward equation ( |2.21|) , for the control history h(t) which 
appears as the argument of Wz- The result is obtained by simply substituting the constraint 
( ggg ) into the inverse Legendre transform for Wz- 



W z [h] = [ * dth T {t)z(t)-r z [z] 
Ju 

= f f dt {h T (t) [(Z(t)) t + (Z(t)C(t)) t ] - < C(t), (d t - L)V{t) >} 

J ti 

= f 1 dt {h T (t)(z(t)) t - < c(t), (d t - L h )v(t) >} 

dth T (t)(Z(t)) t , (2.26) 



The second term in the third line vanishes by equation ( |2.21| ), 

The relation ( 2.25j ) is a compact presentation of the fluctuation-response relations. Let us 



define the response functional of order n for the variable Z(i) as 



where (Z(t))t denotes as before the average with respect to the solution V(t; h) of (2.21). This 



functional is causal, i.e. it vanishes if t < ti for any i = l,...n, because the average (Z(t))t 
cannot have any dependence upon h(i') for t! > t. The nth-order response function is taken to 
be the value at h = 0: 

Ri;ii-i n (t;h, ■■■,t n ) ■= l?i ; i 1 ...j n [t;£i,...,t n ;h]| h= Q . (2.28) 
We now state the fluctuation-response relations: 

Proposition: The (n+l)st-order cumulant Cj 1 ---i n+1 (ti, t n+ \) is determined for each integer 
n > by 

n+l 

Cii-i„+i(*l, ••■>*n+l) = ^2 ^i k ;ii-V k -in+S tk,tl ' — >*n+l)> (2.29) 

fc=l 

where the hat " ^ " denotes omission of the corresponding expression. 

Observe that only the one term with t^ = max p t p is actually non-zero in the sum, by causality 
of the response function. 
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The proof of the relations is very simple. We recall the formula fl2,5Q for the cumulant in 
terms of the generating functional Wz and the corresponding definition Q2.27 ) of the response 



functionals. Taking n + 1 functional derivatives of Wz in ( |2.25 ), one obtains 

S n+1 W z [h] "+ 1 
Sh il (t l )---Sh in+1 (t n+1 ) - ^ ^;i 1 -^-W 1 ^'* 1 '-' tfe '-'*" +1 ' nj 

+ Z) / ^/»i(*)i2j i i 1 ..^ 1 [<;*i ) -,*n+i;h]. (2.30) 
This formula is easy to prove by induction upon n. Setting then h = 0, the last integral term 



vanishes and one obtains (2.29). 



It is important to point out that the operator Lh in (2.21) has a quite simple and explicit 



dependence upon the control field h(t), given in ( p. 15 ). While simple, the coupling does depend 



upon the solution V{t) itself, through the average (Z{t)) t - Hence, the equation ( 2.21 ) for 
V(t) is actually quadratically nonlinear, unlike the original master equation. Nevertheless, this 
nonlinearity is exactly that required to preserve the normalization of the solution. The equation 
may be integrated forward in time to obtain simultaneously V(t) and (Z(t))t as functionals of 
the control h. The response functions can then be determined by differentiating the results. 

3 Moment- Closure Approximations 

(3.1) Variational formulation of moment- closure 

The results of the previous section provide a general approach to computation of the multi-time 
cumulants. However, it is obvious that the required integration of the modified master equa- 
tion ( p. 21 ) will be possible only for the simplest of models, with a few degrees of freedom. For 



spatially-extended systems with many degrees of freedom, this integration is totally intractable. 
The practical employment of the FRR's then depends upon making moment-closure approxi- 
mations. We shall review here the variational formulation of moment-closure approximation, 
following essentially the treatment in However, we shall also introduce some important 
simplifications, that we comment upon as we proceed. 
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Moment-closure approximations to the generating functional Tz[z] of irreducible multi-time 
correlations of Z(t) are obtained by means of the characterization of that functional through the 
constrained variation in ( |2.8| ). Rather than varying over all A £ L°°,V € L 1 , one varies only 
over finitely parametrized trial functions. The trial functions are constructed from the usual 
elements of a moment-closure: a set of moment functions Mi(x,t), i = 1,...,R and a PDF 
Ansatz "P(x, t; fi), which is conveniently parametrized by the mean values which it attributes to 
the moment-functions, fi := J dx V(x, t; /x)M(x, t). The left trial function is then taken to be 

R 

B(x,t;o) :=53oiMi(x,t). (3.1) 

i=i 



Following the discussion in section 2.1, we have chosen the left trial state in the form ( 2.11| ), 



to incorporate automatically the overlap constraint fl2.9|) . The histories a(t),fj,(t) are the 
parameters to be varied over. Substituting the trial forms, one obtains the reduced action 

T[a,ti] = I* dt a T (t)[U(t)-V(n(t),t)} (3.2) 
Ju 

with 

V(/x,*):=<($ + £*)M(t)) M . (3.3) 

Of course, (■)// denotes average with respect to the PDF Ansatz. An unconstrained variation 
of ( p.2j ) recovers the standard moment-closure equation: = V(/i, t). 

For the calculation of the action, however, there is the additional expectation constraint 
( 2.10| ). In terms of the trial parameters, it becomes 



z(t) = C(/x(*),t) + C z (ti(t), t)a(t). (3.4) 

Here, 

CQM) :=<Z(t)>^ (3.5) 
is the Z-expectation within the PDF Ansatz and 

Cs (a*, t) := (Z(t)M T (t)) M - C(A*, t)/^ T (3.6) 
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is the corresponding ZM-covariance matrix. It is remarkable that t), CzQu, t) are the only 
inputs of the PDF Ansatz actually required for the calculation. When the constraint ( |3.4[ ) 
is incorporated into the action functional Q3,2| ) by means of a Lagrange multiplier h(i), the 
resulting Euler-Lagrange equations are 

ii = V(//,i) + C5(/M)h(i) 

:= V z (ji,h,t). (3.7) 

and 

a + (j£f) T (J*, h,t) a +^Y (fi, t)h(t) = 0. (3.8) 

These are solved subject to an initial condition /x(tj) = /x and a final condition oc(tf) = 0. 
When the solutions of the integrations are substituted into (|3.2| ), there results a Rayleigh-Ritz 
approximation Tz[z] to the effective action of Z(i). The value z(t) of the argument is that given 



by the constraint equation (3.4) for the given value of the control field h(t). 

The above derivation of the moment-closure approximation Tz[z] is equivalent to that in 
|Q] but differs in some important details. The trial states employed in y] contained each an 
additional parameter, /iq and ao, with JI = /iq(1, fi),a = (ao,a). Thus, the trial states 
employed there may be written, in our present notations, as 

7>(x,t;/Z) :=mP(xi,t;n) (3.9) 

and 

R 

A(x,t;a) :=^aiMi(x,t). (3.10) 
%=o 

Thus, fiQ was an arbitrary normalization factor and ao was the coefficient of the constant 
moment-function Mo(x, t) = 1. With this pair of trial functions, the overlap constraint ( p.9[ ) 
was no longer automatically enforced and needed to be incorporated via a Lagrange multiplier 
X(t). The resulting Euler-Lagrange equations of the constrained variation for ~p(t),a(t) then 
involved both multipliers h(t) and X(t). See the equations (3.93)-(3.95) in H. Nevertheless, 
those equations are equivalent to (|3.7|)-(|3T8|) above. We shall not give all the details here, but 
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leave it as a relatively simple exercise for the reader to check. We only point out that one may 
take always fiQ = 1, without any loss of generality, by absorbing that factor into the coefficients 
a of the left trial function. It then follows from the 0-component of the equation (3.93) in |3|] 
that the Lagrange multiplier for the overlap constraint is given explicitly as 

A(t) = Fo(/x, h, t) = h T (t)C(/x, t). (3.11) 

If one uses this result and also uses the constraint equation (3.95) in to eliminate the variable 
cto from the equations, then one derives from (3.93)-(3.95) in Q identically the same equations 
as Q-(]3j) above. 

Despite their equivalence to the variational equations in Q, the new form in ( [3.7D -(p\8|) 
above are far more convenient. Because of the presence of the multiplier X(t) in both the 
forward and backward equations (3.93)-(3.94) in Q, those equations posed — apparently — a true 
initial-final value problem. It was proposed in Q to solve that boundary- value problem in time 
with a relaxation method. However, we now see that the forward equation ([Tt]) is completely 
uncoupled from the backward equation. It may be integrated forward in time, storing the 
solution fi(t) for subsequent input into the backward equation ( [3.8D for a(t). Efficient numerical 
algorithms for doing so and then integrating the results to calculate the approximate action 
have been developed by us and will be discussed in another work. 
(3.2) Fluctuation-response relations in closures 

A Rayleigh-Ritz approximation Wz[h] to the cumulant- generating functional may be introduced 
by the formal relation 

W z [h}+f z [z] =<h,z> (3.12) 
It may be easily checked that this definition is equivalent to the formula 

W z [h] = f tf dt h T (t)t(ti(t),t) (3.13) 

Jti 

in which /x(i) is the solution of the forward equation ( [3.7D for the control history h(t). The 
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derivation is the exact analogue of ( |2.26| ). Indeed, it follows that 

rt f 



W z [h] = [ f dt h T (t)z(t)-f z [z] 

hi 

h dt {h T (t) [C(/x(t),t) + C z (fx(t),t)cx(t)} - a T (t)[A(t) - V(/x(t),t)]} 

i 

*' dt {h T (t)£(Mt),t) - <* T (t)[Kt) ~ Vz(M(t), h(t), t)]} 



= [''dt h T (t)C(v(t),t), (3.14) 

where ( |3.7|) was used to eliminate the second term of the third line. 

It is a very attractive feature of the Rayleigh-Ritz approximation scheme that the resulting 
functionals T z [z], W^[h] remain formal Legendre transforms of each other. That is, 

ST 7 5W 7 

This follows from the discussion in 0, where it was derived by a reduction from the under- 
lying variational formulation of the master equation. It is worthwhile to record here another 
derivation, which is instead based directly upon the constrained moment equations ( |3.7D -(|3l|). 
Among other things, this new proof carries over usefully to discrete approximations of the 
moment equations employed in numerical computations that do not follow directly from an 
underlying microscopic theory. 

The derivation begins by functionally differentiating ( |3.13| ) with respect to h(t): 



Sh(t) 



ds [im) h(s) - (3 - 16) 



The abbreviation := £(fi(t),t) was introduced and in the second line the chain rule was 
employed. Causality was also invoked to reset the lower limit of integration. The response 
matrix S g^us that now appears satisfies an equation obtained by functionally differentiating 



with respect to h(t): 
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with the abbreviation 



A(s) := ^-(fx(s),h(s),s). 



(3.18) 



The equation ( |3.17 ) can be solved with a Greens function given by a time-ordered exponential: 



Sh(t) 



Texp 



A(r) dr 



(3.19) 



The adjoint of this propagator appears in the solution of the backward equation (|3.Sj). That 
equation can be written as 



d(t) + A T (i)a(i) + (J£ ) T (/i, t)h(t) = 



(3.20) 



and its solution is 



a(t) 



■f — 



ds Texp 



A (r) 



— ) (^..s)h(.S-) 



(3.21) 



Here Texp [■] denotes the anti-time-ordered exponential. If the solution ( 3.19 ) for the response 
matrix is substituted into the second term of fl3.16|) , and then ( |3.21 ) is employed, it follows that 

5W Z 



5h(t) 



[h] = C(t) + C z (t) a (t) = z(t), 



(3.22) 



using (|3.4j) . Thus, the second relation in ( |3.15 ) is proved. Of course, the dual first relation is 
obtained by functionally differentiating ( |3.12[) with respect to z(t) and using ( |3.22| ). 

The Rayleigh-Ritz approximate generating functional W^^fh] given in ( 3.1 3| ) retains most 
of the remarkable features of the exact generating functional W^[h]. In particular, its value 
may be obtained by integrating just the forward moment-equation ( |3.7D for the selected control 
field h(i), and then substituting the solution /j,(t; h) into ( 3.1 3| ) . Fluctuation-response relations 
follow for the approximate cumulants in the same way as before. Just as before, one may define 
the approximate response function of order n for the variable Z(t) as 



^i;ivi n ^1) •••> ^n) 



5h h (ti) • • • 5h in (t n ) 



h=0 



(3.23) 
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This functional is also causal. By the same argument as before, one obtains the fluctuation- 
response relations for the Rayleigh-Ritz approximate cumulants generated from W^[h]: 



n+l 

C il ... in+1 (t 1 ,...,t n+ l) = R i k ;i x . 

fc=l 



{tk'i *1> - j tki •■■! in+l)- 



(3.24) 



■*fc'"*n+l 

This FRR within moment-closures is a practical way to compute multi-time correlations nu- 
merically, as we shall see below. 
(3.3) Physical Interpretation of the closure FRR 

The results are easiest to interpret physically in the case n = 1. In that case, the FRR deals 
with the 2nd-order cumulant of Z(t) or the 2-time covariance matrix C(t, to) := (5Z(t)5Z T (to)) 
(with 5Z(t) := Z(t) - (Z(t))). The FRR here states that 

C(t, t ) = R(t, ^o) + [R(to, t)] T , (3.25) 



with R(Mo) :-- 



6h(t ) 



h=0 

C(Mo) 



, the response matrix. Thus, for t > to, 
5C(t) 



h=0 



Sh(t ) 

(a t) ML 
(M^)'Texp 



h=0 

t 



to 



A*(s) cfe 



CjGMo) 



(3.26) 



using ( pTT9|) . We set A*(t) := A(i)| h=0 . Recall also that C^(/J,i ) = (<5M(i )<5Z 1 (t )). 

We see that the same result can be obtained by making two physically-motivated approxi- 
mations. The first is the slaving hypothesis: that fluctuations of the variable Z(t) are instanta- 
neously slaved to those of the moment-variables M(i), or 5Z(t) ~ ( ^ J (fi, t)5M(t). Thus, 



(5Z(t)5Z T (t )) « ( |£) (/x,t)(5M(t)5Z T (t )). 



(3.27) 



The second approximation is the regression hypothesis: that fluctuations of the moment- 
variables M(t) decay on average according to the linearized closure equation, <5M(i) pa A*(i)<5M(i). 
Thus, 



{5M(t)5Z T (t )) w Texp 



to 



A*(s) (is 



(£M(t )5Z T (* )). 



(3.28) 



Together, (ggg ) and ( ^28|) lead directly back to the result <^2^) . 
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The special case of the FRR in (3.25) with n = 1 and with Z(t) taken to be the moment- 
variable M(t) itself, was previously derived in ||. It was already pointed out there that the 
physical interpretation of the approximate FRR was provided by the regression hypothesis. 
That result has now been generalized to the case where Z(i) ^ M(i) and the new physical 
principle in the Rayleigh-Ritz approximation is the slaving hypothesis. 

Of course, it is also of interest to consider the physical meaning of the approximations 
involved for n > 1. For the next case, n = 2, the object of interest is the 3-order cumulant 



Cijk(t2,ti,t ) := {5Zi(t2)5Zj(ti)SZk(to)). 



(3.29) 



We consider, without loss of generality, the case t<i > t\ > to- Using the FRR and the chain 
rule twice, we obtain 



Cijk{t2,h,to) 



h=0 



8hj(ti)5h k (to) 

d% ( } 



<W(*2) 

h=0 Sh k (t Q ) 



h=0 



+ d lM [ ^ t2) 5hj{t x )5h k {to) 



(3.30) 



h=0 



We see that the slaving principle holds in a generalized sense. Now higher order derivative 
terms in the Taylor expansion of C(/x, t) appear beyond the leading one. 

In order to focus on the dynamical aspects of the Rayleigh-Ritz approximation for n = 2, 
let us consider now just the special case Z(t) = M(t), so that £(fj,,t) = [i. We shall show that 
the FRR for n = 2 (and Z(t) = M(t)) implies that 



Cijk(h, h,t ) 



*2 Q 2 V 

dt Ei p (t2,t)— — (fl, t)C r j(t, ti)C qk (t, to) 



+ Ei p (t2,ti, 



dC~ 



jp 



h)C q k(ti,to). 



(3.31) 



We have introduced the propagator of the linearized closure dynamics: 



E(M') := Texp 



A*(s) ds 



(3.32) 
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The result ( |3.31| ) is obtained by taking the second functional derivative with respect to h(t\] 
of the lst-order response functional 



R[t2,*o;h] = Texp 



A(s) ds 



t 



CQMo) 



(3.33) 



using the simple identity (a continuous "product rule" of functional differentiation) 



5 



5hj(t\ 



-Texp 



t2 



A(t) dt 



/ti r rt2 

dt T exp J A(s) ds 



5A(t) 



T exp 



A(s) ds 



to 



(3.34) 



and then setting h = 0. 

Some physical insight into the result ( |3.31 ) of the closure FRR can be obtained by rederiving 
the result in a more heuristic way. Let us make a nonlinear regression hypothesis: that the 
fluctuations evolve, in general, according to the full closure dynamics. By a Taylor expansion 
to quadratic order, one then obtains 



SMi(t) = A^jtySMjit) + \ d d ^ (iM,t)5M q (t)5M r (t) + O (<5M 3 ) . 

This equation can be rewritten in integral form as 

1 r t2 .. „ ,. d 2 V r 



(3.35) 



6Mi(t 2 ) = Eipih^SMpih) + - f 2 dt E ip (t 2 ,t) 



dfiqdflr 



(v,t)5M g (t)5M r (t) + O (5M 3 ) . 

(3.36) 



Multiplying by 6Mj(ti)6Mf.(to) and averaging then yields the formula 
(SM^SMj^SMkito)) ; 



I f Z dt E ip (t2,t)^-^( t i,t){5M q (t)5M r (t)5M j (t 1 )6M k (t )) 

l Jti O^gOflr 



+E ip (t 2 ,t 1 )(5M p (t 1 )5M j (t 1 )6M k (t )). 



(3.37) 



Next, one can approximate the remaining correlations by discarding all cumulants of higher 
order than third (and disconnected terms single-time in t). Then, 

(SM^SMrityMjitJSMkito)) « C qj {t,t x )C rk (t,t Q ) + (q «-> r), (3.38) 

which relation, substituted into the first term of (3.37) yields precisely the first term of ( 3.31| ). 
Likewise, 

{SMp^SMjitJSMkita)) » -Cpi^C^^Y^t^C^txM), (3-39) 
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where Ti mq (fi,ti) is the instantaneous 3rd-order irreducible correlation in the PDF Ansatz. 
Substitution into the second term of ( |3.37| ) yields the corresponding term of ( 3.31| ) if we assume 
that 

dC 

-rr^Ox, *i) = -C p i{t x )C jm (tx)Y lmq {ti) (3-40) 

instantaneously at time t\. 

Thus, the two terms in ( |3 .31 ) have quite different physical interpretations. The first integral 



term represents triple correlations dynamically generated at intermediate times ti > t > t\ from 
the correlations propagating in from times t\ and to, which subsequently then relax to time t?. 

Q2y 

according to the linearized closure dynamics. The second-derivative dfi can be interpreted 
as an "interaction vertex", due to the nonlinear terms in the closure equation, by means of 
which the fluctuations interact. The second term, on the other hand, is not produced by any 
interaction of the fluctuations. It represents a triple correlation present instantaneously in the 
closure PDF Ansatz which is simply propagated in time by the linearized dynamics. Both of 
these terms, as well as the terms for cases n > 2, can be expressed in terms of suitable diagrams. 
These involve propagators from the linearized dynamics and vertices both from the nonlinear 
terms in the dynamics and the higher-order correlations in the instantaneous PDF Ansatz. 
See the Appendix. Of course, these are not "bare" Feynman diagrams, for the Rayleigh-Ritz 
approximation is highly non-perturbative and the propagators and vertices represent "dressed" 
objects resulting from statistical closure. 

The single-time relation ( |3.40| ) that we derived heuristically is, in fact, a necessary condition 
for consistency of the Rayleigh-Ritz approximation to the triple correlation. Only if it is true 
will the 3-time correlation in ( |3.31| ) coincide along the diagonal ti = t\ = to = t with the 
value Cijk(n,t) calculated from the single-time PDF Ansatz V(x; fi,t) which is input into the 
Rayleigh-Ritz calculation. Indeed, if we assume that ( |3.40| ) holds, then the second term of 



(3.31) can be rewritten as 

Cijl(t2,h,to) = —Cii{t2,ti)Cj m {ti)T lmq {ti)C q k{ti,to), (3-41) 
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from which it is manifest that Cijk(t, t, t) = Cjjfe(/x, t). Relations such as ( |3.40| ) are also very im- 
portant in other contexts within the Rayleigh-Ritz method. For example, a relation equivalent 
to ( fUOD , or 

r ijk (p,t) = —^(n,t), (3.42) 

was employed in fla] to prove an .ff-theorem at quadratic order. 

However, the relations ( |3.40| ),( |3.42 ) are not automatically true for an arbitrary PDF Ansatz. 
They might be taken as definitions of the triple-correlations within the Rayleigh-Ritz approx- 
imation, which, we should remember, has available to it directly from the single-time PDF 
Ansatz only the mean ((fi,t) = /x and the covariance Cz(fJ>, t) = C(/x,i). However, the triple 
correlators Cijk(n,t) and Tijk(fi,t) are symmetric in their indices i,j,k, while the definitions 
through ( |3.40| ) and ( p. 42 ) need not have such symmetry. Thus, such a definition may not be 
consistent. Fortunately, it is possible to construct the closure to ensure that (|3l0| ), (|3l2|) 
hold, by employing an exponential PDF Ansatz, such as those previously developed for Boltz- 
mann kinetic equations in transport theory ||. See also ||. Within such a closure scheme the 
single-time irreducible correlators are all obtained from a generating function: 

T il ... in ( t i,t) = (aM) (3.43) 

In fact, H(fi,t) is just the relative entropy. Because of ( |3.43 ), the consistency condition ( 3.42| ), 
as well as all higher-order ones, can be automatically ensured by constructing the closure via 
an exponential PDF Ansatz. Thus, this closure methodology has a special relation with the 
Rayleigh-Ritz approximation scheme. This will be discussed in detail elsewhere Wn. 



4 The FRR in Numerical Computation 

(4-1) Numerical Differentiation 

We have studied some of the properties of a PDF Ansatz required for a physically accurate 
and mathematically consistent approximation of time-correlations via the FRR's. Another im- 
portant issue is the feasibility and accuracy of the FRR's for use in numerical computations. 
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Except in special circumstances, it will not be possible to employ the FRR's analytically and 
numerical solution on the computer will be required. We have seen that the FRR's give the 
time-correlations by (functionally) differentiating the solutions of the modified closure equation 
with respect to the control field h. The numerical problem is to compute the required 
derivatives. It is well-known that finite-difference approximations of derivatives are inherently 
numerically unstable, because the decrease in differentiation step Ah needed to reduce trunca- 
tion error must cause the round-off error in finite-precision arithmetic to grow. If the FRR's are 
to be a useful computational tool, better numerical differentiation methods must be devised. 

Fortunately, this problem has been encountered and solved in the context of other dynamical 
problems. One of the main application areas is sensitivity analysis, in which the sensitivity of the 
solution of a dynamical equation to changes in the initial data or to parameters in the equation is 
required []|]. By "sensitivity" we mean just the Jacobian derivative matrix of the solution vector 
with respect to the parameter vector (or higher-order derivatives). Our problem is exactly of 
this form, in which the "sensitivities" required are those of the solution of the modified closure 
equation with respect to the added control field h. The numerical techniques which yield 
accurate derivatives in sensitivity analysis depend upon solving additional dynamical equations 
for the derivatives themselves. There are two general techniques for doing so, depending upon 
the time-order of propagating derivatives: the "forward mode" and the "reverse mode" . These 
two techniques have already been illustrated in the context of our earlier discussion. The 
equation ( |3.17 ) for the response matrix is equivalent to 



Sh(t ) Sh(t 
integrated forward in time with initial data 

5/i(t) 



6m A(t)M (4.1) 



Sh(t ) 



= Ck(p,to)- (4-2) 

t=t 



Substituting the result into ( 3.16 ) gives the derivative sh(t Z ) M • This illustrates the forward 



mode. Alternatively, one may compute the same derivative by integrating the adjoint equation 
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720| ) for ot(t) backward in time and then substituting the result into the formula 

^[h]=C(t) + C z (t)a(t). (4.3) 

This illustrates the reverse mode. Such equations may be developed for arbitrary derivatives 
and, implemented numerically, they yield accurate and stable approximations. 

Not only are these approaches numerically efficient but they can also be largely automated. 
Software for "automatic differentiation" is now becoming widely available: see Q. Such tools 
directly generate from source code for numerical computation of the solution of the dynamical 
equation a corresponding code for the computation of its derivative. There is no need to 
compute required input derivatives, such as Ay = by hand. Furthermore, it easy to 
compute "sensitivities" with respect to new perturbations, such as those corresponding to a 
new class of variables Z(t) of interest, without requiring extensive recalculations. The method 
has been tested and proved successful in application to real-life codes for PDE's employed in fluid 
dynamics and elsewhere. The availability of such software greatly enhances the attractiveness 
of the FRR's as a computational method to calculate time-correlations. 
(4-2) A Simple Example 

To illustrate the computational use of the FRR, we shall consider a simple closure for the decay 
of homogeneous, isotropic turbulence governed by the Navier-Stokes dynamics. This closure was 
originally employed by Kolmogorov to predict the mean energy decay. It was employed within 



the Rayleigh-Ritz formalism in [10| to predict the 2-time correlation of the energy fluctuations. 
It should be emphasized that this closure omits a physical effect which is very important in 
the decay of energy fluctuations: their relaxation by turbulent diffusion in space. To see such 
effects, one must construct the closure not just for the kinetic energy at a single point (say, 
the origin) but with the kinetic energy at all space points as closure variables. In that case, 
the closure equations contain "eddy viscosity" terms, which are an important linear relaxation 
mechanism for fluctuations. Such improvements have been investigated and tested against 



simulation data in [11]. However, the 1-moment closure is adequate for our purpose, which is to 
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study the utility of the FRR for numerical computations. The main merit of the closure is that 
the Rayleigh-Ritz 2-time correlation is given analytically: see Eq.(4.4) in [1C]. This provides 
an objective basis of comparison for numerical results. We shall only consider here the 2-time 
correlations, i.e. the case n = 1. 

The closure we consider has just one moment function, the kinetic energy K = \v 2 at a 
single point in space. The moment average \x is here denoted E. It obeys the equation 



E(t) = -A m E p (t) 



(4.4) 



in which A m ,p are suitable real constants. See [1C] for details. The single-time covariance 
C{t) := {[5K(t)Y) is given in the closure by 

2 



C(E;t) = ^E 2 , 



(4.5) 



which follows from assuming a Gaussian 1-point velocity distribution. Thus, the perturbed 
closure equation becomes here 

E = V(E) + h{t)C(E) (4.6) 



in which V(E), C(E) are given by (f4.4|),( fI~5D , respectively. Now the FRR states that the 2-time 
correlation C(t,to) := (6K(t)SK(to)) is given in the Rayleigh-Ritz approximation by 

6E(t) 



C(t,*o) 



(4.7) 



h=0 



Sh{t ) 
for t > to- 

The righthand side of ( [4.7D has been calculated by us numerically, in two different ways. 
The first method is based upon the observation that 



SE(t) 



Sh(t ) 



h=0 



dE 
~dh 



(t;h) 



h=0 



(4.8 



where E(t; h) is the solution of the closure equation for the modified initial data 



E(t ;h) :=E + h-C(E ;t ). 



(4.9) 
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The closure equations were numerically integrated with a 4th-order Runge-Kutta scheme with 
time step At = 10 -3 , in double precision arithmetic, but with initial data given by ( |4.9| ) for 
the two small values h = ±10 -6 . The derivative ( (4.8| ) of the solution at later times t was 
then estimated by the symmetric, 2nd-order finite difference approximation to the derivative. 
The second method for calculating the functional derivative in Q4.7D is to solve the analogue 
of equation ( fP| ) with initial data t t = ^^ 0, ^°^ 1 together with the closure equation 

itself. These were both integrated numerically by the same Runge-Kutta code as before. The 
matrix A(i) that appears in fl4.1|) (here, a 1 x 1 matrix) is given analytically by 

A,(t) = - P A m EP-\t) (4.10) 

and it was input directly into the code. Hence, the only errors in the functional derivative 
calculated by this second method, as in the solution of the closure dynamics for the means, are 
the 4th-order truncation errors and the round-off errors. 

We show in Figure 1 the correlation C(t,to) calculated by the FRR, compared with the 
analytical Rayleigh-Ritz solution given in Eq.(4.4) of flQfl . Both the values calculated by the 
finite difference approximation (method 1) and the adjoined equation for the Jacobian (method 
2) are shown. As may be seen, these agree perfectly, both with each other and with the exact 
result. This is not to be unexpected in our example, considering the high-order accuracy of 
our approximations and the double precision arithmetic. In other cases, it is hard to assess a 
priori the accuracy of the finite difference approximation, for which the optimal discretization 
step-size Ah op t may depend upon both space and time in an unknown manner. 

It was shown in || that, in general, the same 2-time correlations provided by the FRR are 
also given by a linear Langevin equation. That model may be only formal, since the covariance 
of its noise term need not be positive. In any case, the numerical use of the Langevin model to 
calculate the 2-time correlations is far less efficient than the numerical use of the FRR. Not only 
must the stochastic equation be integrated for a large enough number of realizations N 1, 
but also in each realization random number generators must be called in each step of the time 
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integration. Furthermore, the individual realizations governed by the Langevin dynamics shall 
be far less smooth in space and time than averages over the ensemble, and thus much smaller 
space and time discretization steps Ax, At will be required. The computational expense of 
using the Langevin model is, thus, far higher than the FRR and is not to be recommended. 
The stochastic equation is only useful for conceptual purposes. The FRR, on the other hand, 
is quite efficient because it takes full advantage of the increased regularity and stability of 
statistically-averaged quantities. It is really a "thermodynamic approach" to calculating the 
time-correlations and not a "statistical-mechanical" method. 



5 Conclusions 

In this paper we have reviewed and simplified the variational approach to statistical dynamics 
proposed in As a main new result, we derived a general fluctuation-response relation 
(FRR) for arbitrary multi-time correlations. We demonstrated that the FRR's are preserved 
in a moment-closure approximation by the Rayleigh-Ritz method. We discussed the physical 
significance of the closure FRR's in terms of various intuitive hypotheses: slaving, regression 
(linear and nonlinear) of fluctuations. We also discussed computationally efficient and accurate 
methods for computing the derivatives required in the FRR's. 

Many interesting problems can be investigated with the present methods. These include 
temporal multiscaling in turbulence fll2| [H^], aging phenomena in glassy relaxation fi"4| , p~5[| , 



transition rate theory in chemical kinetics [16], and Lagrangian statistics of advected scalar 



reactants [|17]]. The FRR should also hold for quantum systems, governed by quantum Liouville 
or master equations. Rayleigh-Ritz methods could provide a tractable means to compute multi- 
time statistics in quantum field theory and in the quantum many-body problem. 
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Appendix: Diagrammatic Rules 

We sketch concisely here the diagrammatic rules which follow from the closure FRR for the 
multi-time cumulants or connected correlation functions. The derivation generalizes that for 
the 3-time cumulant in subsection 3.3. It is advantageous to specify the latest time to be t n 
and the earliest time t$. Then, (n + l)-time correlations for tf = t n > t n _x > • • • > t% > to = ti 
are obtained by succcessive functional differentiation of expression ( gjgD for R[t n ,t ;h] with 



respect to h(ii), h(t n _i), using the "product rule" Q3.34p . We recall from ( |3.18| ) that 

, . , <9V , . , T , , dC . . 

and employ the chain rule to calculate successive derivatives. We also assume that the single- 
time irreducible correlations have the entropy as a generating function, as in ( 3.43| ). 

The terms which result may be associated to graphs. The lines in the graphs terminating 
at times t, t? (internal or external) are given by the covariance function Cij(t,t'). If t = t', then 
Cij(t,t) = Cij(fi,t). The vertices are of two types. For each integer r > 2 there are (r + l)-fold 
vertices of the form 

Wf ir (s) :=T im (8) ^ {H,s) 

and 

d r+1 H 

- r ii-ir+i(*) = ~7j- —(/*.*). 

U h l jl ' ' ' H'jr 

where the latter is just minus the single-time, irreducible (r + l)st-order correlator. The minus 
sign appears because of the fact that C = T 1 and thus 

One may replace the VF-vertices with V^ " ^ r {s) := d ® Vm ^. Qu, s) if the propagator line C^s', s) 
entering the i-node of the VF-vertex is replaced by a linear propagator Ek m (s',s) entering the 
m-node of the V- vertex. 

The following rules apply: 
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(i) The graphs which appear are all tree graphs with the times t n ,t n -i, ...,ti,to terminating 
the external lines. The trees are rooted at time t n and branch up to earlier times, with 
the times non-increasing as one ascends the tree. 

(ii) Each vertex must be linked to at least one of the early external times t„_i, ...,io directly 
by a propagator line. 

(iii) The T-type vertices are all evaluated at an early external time t n -±, ...,to which is deter- 
mined as the latest time reached by any branch starting upward from that vertex and 
passing only through T-type vertices. 

(iv) The W-type vertices (or the F-type) are evaluated at internal times s which are integrated 
over the largest possible subrange of t n > s > to consistent with the rule of non-increasing 
times ascending the tree. 

Because of rule (ii) it is clear that there are only finitely many graphs contributing to each 
(n + l)-time cumulant function, with vertices of at most (n + l)st-order appearing. The finite 
sum of all the contributions from these graphs gives the FRR result for the cumulant function. 
Thus, it is clear that this graphical representation is not a perturbation expansion into Feynman 
diagrams, since the latter contains closed loops and infinitely many terms. The propagators 
and vertices here are all "dressed objects" and the representation is nonperturbative. 

When all the external times are equal, t n = ■ ■ ■ = to = t, then the graphical representation 
simplifies considerably. There are then no V or W-type vertices, because the integration range 
over each internal time s has shrunk to zero. Furthermore, all of the T-type vertices are now 
evaluated at the same time t. In fact, the resulting graphical expansion is just that of the 
well-known representation of the single-time (n + l)st-order cumulant Ci 1 ...j n+1 (t) as a sum over 
tree diagrams with single-time irreducible correlations Ti 1 ...i r+1 (t) as vertices and 2nd-order 
correlators Cij(t) on the internal and external lines. Thus, we obtain a proof for any order 
(n + 1) that, along the diagonal t n = ■ ■ ■ = t = t in time, C il ... in+1 (t, ...,t) = C il ... in+1 (p,,t). 
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Figure 1: Comparison of Two-Time Correlations. FRR Calculations vs. Exact Result. 
Initial energy at to = was set to Eq = 10. The decay was calculated with spectral parameters: 
2, A = 1, and Kolmogorov constant a = 2. Notations follow reference [10]. 
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